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ABSTRACT 


^This  paper  addresses  the  issue  of  load  path  discretization  in  in¬ 
cremental  tangent  modulus  finite  element  analysis.  The  view  is  taken 
that  discretization  can  best  be  accomplished  during  the  course  of  the 
numerical  solution  by  employing  a  constraint  condition  which  restricts 
the  level  of  structural  stiffness  approximation  at  each  load  step.  By 
identifying  a  field  variable  which  strongly  influences  the  stiffness,  a 
suitable  constraint  condition  can  be  selected  in  terms  of  the  nodal  vari¬ 
ables.  With  such  a  constraint,  design  of  a  solution  algorithm  for  deter¬ 
mining  the  step  size  along  with  the  nodal  variables  is  straightforward, 
as  is  demonstrated  in  the  text.  This  variable  load  step  solution  approach 
provides  the  analyst  with  a  simple  yet  efficient  method  for  logically 
controlling  step  size,  without  having  to  resort  to  time  consuming  and 
costly  reanalysis  procedures  to  insure  that  the  numerical  approximation 
is  within  satisfactory  tolerance.  A 


Stepwise  nonlinear  tangent  mod\Tus  formulations  which  employ  average 
stiffness  matrices  for  a  step  are  examined  in  detail.  The  iterative  na¬ 
ture  of  the  solution  method  is  discussed  in  general  terms,  while  specific 
consideration  is  given  to  the  nonhardening  Prandtl-Reuss  elastic-plastic 
problem.  For  this  problem  the  deviatoric  stress  change  which  occurs  dur¬ 
ing  plastic  deformation  is  identified  as  the  crucial  stiffness-governing 
incremental  field  variable.  Test  solutions  are  presented  for  a  constraint 
condition  which  restricts  changes  in  this  variable  to  a  specified  fraction 
of  the  yield  stress.  The  numerical  results  demonstrate  the  propriety  of 
this  particular  condition,  and  the  viability  of  the  variable  load  step 
solution  approach. 
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INTRODUCTION 


A  common  structural  analysis  problem  is  to  determine  the  stress  and  strain 
distribution  throughout  a  structure  at  distinct  stages  of  an  applied  load  history. 
For  incremental  problems,  such  as  the  following  elastic-plastic  flow  theory  ex¬ 
amples,  it  is  necessary  to  trace  the  solution  history  in  discrete  steps  along  the 
load  path  to  the  load  states  of  interest.  In  this  paper  we  are  concerned  with 
accurate  methods  for  discretizing  the  load  path  in  finite  element  tangent  modulus 
analysis . 

A  characteristic  of  incremental  finite  element  formulations  is  that  the 
structural  stiffness  continually  changes  along  the  load  path,  corresponding  to 
the  changes  which  occur  in  those  stress  and  deformation  variables  which  enter  the 
stiffness  definition.  Hence,  for  nontrivial  problems,  finite  load  steps  imply 
that  the  stiffness  matrix  can  never  be  exact  over  the  duration  of  a  step.  The 
solution  at  a  particular  load  will  always  be  dependent  upon  the  prior  load  path 
discretization,  and  the  analyst  must  consider  means  to  measure  and  control  this 
dependency.  It  is  this  specific  form  of  numerical  error  that  concerns  us  here. 

Most  established  finite  element  computer  codes  require  a  priori  specification 
of  the  various  load  steps  of  an  analysis  based  upon  the  best  judgment  of  the  ana¬ 
lyst.  Successive  reanalysis  is  often  undertaken  with  step  size  refinements  accord¬ 
ing  to  the  observed  solution  discretization  dependency.  It  is  clear  that  this 
error  control  strategy  can  be  very  costly  and  inefficient:  it  does  not  directly 
address  the  source  of  the  error,  and  thus  in  general  an  acceptable  solution  can 
be  obtained  only  in  a  haphazard  fashion. 

In  the  present  work  the  view  is  taken  that  discretization  error  can  be  best 
controlled  by  treating  step  size  as  a  variable,  and  determining  it  numerically 
along  with  the  incremental  nodal  variables  by  employing  an  error-regulating  con¬ 
straint  condition  in  conjunction  with  the  usual  stiffness  equilibrium  equations. 
Since  the  approximate  nature  of  the  stiffness  is  the  source  of  the  error,  the  most 
proper  constraint  condition  is  one  which  regulates  the  stiffness  approximation 
level.  In  the  approach  taken  a  particular  form  of  constraint  is  used  which  ful¬ 
fills  this  function  by  regulating  changes  in  a  field  variable  which  is  known  to 
strongly  influence  the  stiffness  definition.  Although  we  have  limited  our  devel¬ 
opment  efforts  to  the  design  and  testing  of  a  solution  algorithm  for  stress  anal¬ 
ysis  of  Prandtl-Reuss  elastic-plastic  materials,  the  variable  load  step  approach 
is  discussed  in  general  terms  since  it  appears  to  be  applicable  to  other  important 
nonlinear  problems  which  are  routinely  treated  by  the  tangent  modulus  method,  such 
as  large  deflection  analysis. 

Since  we  employ  a  constraint  involving  an  incremental  field  variable,  regard¬ 
less  of  the  form  that  the  constraint  takes  or  the  specific  variable  involved,  in 
the  assumed  displacement  finite  element  method  it  is  in  essence  always  a  condition 
on  the  nodal  displacement  change  vector  AU^ ,  for  step  number  d^,  and  it  can  be 
represented  in  general  form  as 

g  (AUi)  *  0.  (1) 

The  actual  constraint  condition  chosen  for  the  nonhardening  Prandtl-Reuss  elastic- 
plastic  examples  limits  the  yield  surface  deviatoric  stress  change  to  a  given 
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fraction  of  the  yield  stress.  For  a  large  deflection  problem  it  might  be  proper 
to  constrain  a  certain  displacement  gradient  distribution.  Regardless,  with  the 
selection  of  the  constraint  condition  (1),  load  path  discretization  can  become  an 
integral  part  of  the  numerical  procedure.  A  solution  algorithm  has  been  developed 
for  the  class  of  tangent  modulus  formulations  which  employ  average  stiffnesses  de¬ 
fined  in  terms  of  the  undetermined  solution,  and  as  a  result  are  stepwise  nonlinear. 
The  approach  can  be  most  readily  appreciated  by  first  considering  an  elementary 
stepwise  linear  formulation  which  has  the  stiffness  matrix  K_  established  on  the 
basis  of  the  initial  state  alone. 

For  proportional  loading  to  final  load  vector  P  the  stepwise  linear  problem 
at  step  i  takes  the  form 

K  AUt  =  XtP  (2) 

where  represents  the  load  step,  with  Aj  being  the  scalar  defining  the  step 

size.  We  wish  to  select  A^  so  that  AUj  satisfies  both  (1)  and  (2).  This  can  be 

readily  accomplished  by  taking  an  arbitrary  value  for  Aj,  solving  (2)  for  AU; , 

and  then  scaling  AU,  and  A^  by  that  factor  needed  to  satisfy  (1).  This  advances 

the  solution  to  the  next  step  and  the  process  is  repeated  until  the  final  load  is 
reached.  In  fact,  some  of  the  early  elastic-plastic  tangent  modulus  formulations1*2 
took  exactly  this  form.  Their  constraint  condition  in  effect  restricted  stress  in¬ 
creases  to  control  the  change  of  the  elastic-plastic  boundary.  We  will  discuss  the 
inadequacies  of  linear  formulations  along  with  corrective  techniques  which  lead  to 
stiffness  averaging  over  a  step. 

We  examine  then  those  formulations  which  employ  a  stiffness  which  is  depen¬ 
dent  upon  the  yet  to  be  determined  vector  AU^ .  The  equilibrium  equation  (2)  now 
becomes  nonlinear  and  takes  the  form 

K  (A Uj)  AUi  =  AiP.  (3) 

Most  of  the  established  methods  use  some  iterative  scheme  for  solving  (3)  while 
holding  Aj  fixed.  The  AU ;  is  obtained  by  successively  solving  for  trials  Alij 
using  a  stiffness  matrix  based  upon  some  previous  estimate,  K  (Ay]'1)  .  The  vari¬ 
able  load  step  algorithm  suggested  here  scales  the  solution  and  Aj  after  each  it¬ 
eration  until  All]  is  found  which  satisfies  both  (1)  and  (3). 

Elastic-plastic  test  problems  were  solved  using  the  stepwise  nonlinear  vari¬ 
able  load  step  algorithm.  The  convergence  properties  of  the  iterative  process 
were  found  to  be  essentially  the  same  as  the  fixed  load  algorithm.  Whereas  it  is 
always  true  that  convergence  rate  slows  as  stiffness  approximation  level  increases, 
an  important  benefit  of  the  variable  load  step  approach  is  that  convergence  rate 
can  be  regulated  <t  priori  by  adjusting  the  severity  of  the  constraint  condition. 

The  success  of  the  approach  depends  entirely  on  the  identification  of  that  variable 
which  governs  the  stiffness  approximation.  The  numerical  results  show  that  a 
proper  variable  and  constraint  condition  has  been  chosen  for  the  nonhardening 
Prandt 1 -Rcuss  problem. 

1.  POPf  ,  <;.  (J.  The  t /’filiation  of  the  Matrix  Displacement  Method  in  Plane  Plastic-Plastic  IVoMerns  in  Proc.  font'.  Matrix  Methods 
in  Structural  Mechanics,  Wright-Patterson  Air  force  Base,  Al  TDL-TR-66-80,  1%5,  p.  635-654. 

2.  YAMAOA,  Y„  YOSHIMl'RA,  N.,  and  SAKl'RAl,  T.  haslic  Stress-Strain  Matrix  anti  its  Application  for  the  Solution  of  Plastic- 
Plastic  Problems  hv  the  Unite  Plement  Method.  Int.  J.  Mcch.  Sei.,  v.  10,  1968,  p.  343-354. 


VARIABLE  LOAD  STEP  SOLUTION  ALGORITHM 

As  we  have  discussed,  once  we  have  decided  on  a  constraint  condition  (1),  a 
procedure  can  be  developed  whereby  the  given  load  path  is  discretized  as  part  of 
the  numerical  solution.  In  this  section  we  describe  the  algorithm  that  was  de¬ 
veloped  for  stepwise  nonlinear  tangent  modulus  formulations  having  an  equilibrium 
equation  of  the  form  (3).  We  concentrated  on  the  discretization  of  a  particular 
proportional  load  segment  £  of  the  entire  load  path.  The  solution  history  cor¬ 
responding  to  any  previous  loading  is  assumed  known.  The  algorithm  is  an  exten¬ 
sion  of  the  iteration  solution  process  that  was  used  by  Marcal  and  King3  for  fixed 
load  steps. 

The  Marcal  and  King  procedure  involves  at  each  cycle  of  iteration  the  forma¬ 
tion  of  a  trial  structural  stiffness  matrix  followed  by  the  solution  of  the  re¬ 
sulting  linear  matrix  equation.  The  stiffness  trial  at  the  j-th  iteration  cycle 
is  formed,  consistent  with  the  averaging  techniques  of  the  formulation,  according 
to  an  estimated  departure  Auj"'  from  the  current  structural  state,  where  AyH  is 
the  solution  of  the  previous  cycle.  Hence,  at  cycle  j  the  governing  equilibrium 
equation  takes  the  form 

K  (Au{_1)  Auj  =  A.P.  (4) 

The  first  cycle  requires  a  guess  AU°  for  definition  of  the  stiffness  matrix. 
Strictly  speaking,  cycling  must  continue  until  successive  trial  solutions  are 
found  to  be  identical,  so  that 

K  (Auj)  Auj  =  A.P  (5) 

and  therefore  AU]  is  then  the  solution  AU^  for  the  step. 

In  our  variable  load  step  approach  we  adjust  Aj  during  the  course  of  solution 
to  find  that  AU|  which  satisfies  both  the  equilibrium  equation  (3)  and  the  condi¬ 
tion  (1).  This  iterative  process  starts  with  an  estimate  of  the  load  step.  A®,  as 
well  as  with  the  guess  AU°.  For  the  first  cycle  the  matrix  equation  for  AU-  takes 
the  form 

K  (AU?)  AU*  =  A°P.  (6) 

In  general  AU3  will  not  satisfy  the  constraint  condition,  although  there  al¬ 
ways  is  a  scalar  multiple  of  this  vector  which  will.  The  operations  required  to 
determine  the  appropriate  scale  factor  depends  upon  the  nature  of  the  constraint. 
Regardless,  the  scale  factor  is  found  and  the  correspondingly  scaled  displacement 
solution  is  used  as  the  trial  vector  for  the  next  cycle  of  iteration.  The  next 
step  size  trial  follows  from  interpreting  the  scale  factor  as  equal  to  A^/A?,  as 
suggested  by  the  linear  nature  of  (6).  In  general  terms,  the  problem  after  (6)  is 
solved  is  to  determine  a|  which  satisfies 

g  (AU*-A*/A°)  =  0.  (7) 

3.  MARCAL,  P.  V.,  and  KING,  I.  P.  Elastic- Plastic  Analysis  of  Two-Dimensional  Stress  Systems  by  the  Finite  Element  Method. 

Int.  J.  Mech.  Sci.,  v.  9,  1967,  p.  143-155. 
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The  above  operations  for  the  first  cycle  of  iteration  sets  the  pattern  for  subse¬ 
quent  cycles.  At  cycle  j  a  stiffness  matrix  is  formed  according  to  the  estimated 
displacement  AU]'1  •X|'1  /xj'~  .  and  a  new  displacement  AU]  is  determined  from 

K  AU^  =  A  ?  ~  1 P .  (8) 

—  — 1  i  l  — i  l  — 

Once  AU^  is  obtained,  A^  follows  from 

— l  l 

g  (AJjj  *  A ]/a|”  1 )  =  0.  (9) 

We  are  of  course  seeking  in  this  iterative  process  the  step  size  A^  and  the 
associated  vector  AU; .  Our  definition  of  absolute  convergence  must  be  expanded  to 
include  the  step  size  parameter.  Fortunately  we  are  able  to  argue  the  concurrent 
convergence  of  both  AU^  and  Aj.  Thus  the  rate  of  convergence  of  the  iteration 
scheme  can  be  conveniently  monitored  by  a  test  on  the  eye le-to-eycle  change  in  A]. 
The  convergence  test  implemented  in  this  work  states  that  iteration  terminates 
when  the  relative  change  in  A]  in  two  successive  cycles  falls  below  a  given 
tolerance  5. 

We  have  observed  that  the  stress  solution  in  our  elastic-plastic  problems  is 
sensitive  to  the  choice  of  5.  While  this  is  not  surprising,  it  forcefully  demon¬ 
strates  that  the  level  of  convergence  can  govern  the  viability  of  the  numerical 
formulation.  Again,  the  issue  is  most  conveniently  discussed  in  terms  of  the 
fixed  load  step  solution  algorithm.  We  consider  the  ramifications  of  prematurely 
terminating  iteration  at  a  cycle  j,  taking  AU]  as  the  solution  for  the  step.  Re¬ 
ferring  to  (4),  AU]  can  be  interpreted  as  the  equilibrium  solution  for  a  structure 
with  stiffness  (AU]'1)  under  load  AjP.  For  our  elastic-plastic  problems  we  must 
choose  one  of  two  options  for  the  calculation  of  the  resulting  stress  changes  for 
the  step.  The  determination  of  stress  involves  the  use  of  average  constitutive 
relations  based  either  on  AU]'1  or  AU] .  Using  the  first  alternative  for  the  cal¬ 
culation  will  result  in  a  violation  of  the  yield  criterion  and  hardening  law  when 
strain  hardening  occurs.  The  discussion  in  the  next  section  provides  the  back¬ 
ground  for  this  statement.  The  second  option  will  provide  a  stress  solution  which 
satisfies  these  basic  theoretical  constraints  but  is  not  in  equilibrium  with  the 
tractions  represented  by  A|P_.  We  are  unable  to  provide  a  quantitative  relation¬ 
ship  between  overall  load  imbalance  and  our  convergence  tolerance  parameter  6. 
However,  this  is  an  area  which  deserves  careful  consideration  for  both  the  vari¬ 
able  and  fixed  load  step  methods. 

To  complete  our  discussion  of  the  variable  load  step  procedure  we  consider 
some  additional  restrictions  that  should  be  placed  on  the  allowable  magnitude  of 
Aj.  When  a  definite  total  load  vector  I P  is  specified  there  is,  of  course,  the 
need  to  restrict  Aj  <  1.  Furthermore,  EA^  over  all  steps  must  equal  unity.  The 
final  step  to  reach  the  total  load  P  will  usually  be  smaller  than  that  allowed  by 
our  constraint  condition.  For  this  case  the  algorithm  reverts  to  the  standard 
fixed  load  procedure.  When  the  load  vector  P  is  indefinite  in  the  sense  that  the 
final  magnitude  of  its  components  are  not  specified,  then  there  is  no  basis  for 
restricting  the  values  of  Aj.  This  latter  case  applies  to  the  test  problems  con¬ 
sidered  below.  There  the  vector  P  serves  merely  to  specify  load  direction  and  the 
magnitude  increases  step  by  step  without  restriction  until  limit  load  is  detected. 
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CONSTRAINT  CONDITION  FOR  ELASTIC-PLASTIC  PROBLEMS 

The  variable  load  step  solution  method  was  implemented  for  stress  analysis 
of  Prandtl-Reuss  materials.  With  this  material  model  yielding  is  governed  by  the 
Mises  yield  criterion  and  plastic  flow  occurs  according  to  the  yield  surface  nor¬ 
mality  rule.  The  incremental  tangent  modulus  formulation  that  was  employed  has 
been  described  by  Rice  and  Tracey.4  It  is  useful  to  begin  the  discussion  by  con¬ 
sidering  the  general  elastic-plastic  problem  in  the  context  of  various  alternate 
tangent  modulus  formulations.  With  this  as  background,  the  motivation  for  our 
choice  of  constraint  condition  can  be  straightforwardly  explained. 

In  the  elastic-plastic  problem  the  structural  stiffness  changes  as  the 
elastic-plastic  boundary  moves,  and  as  the  direction  of  plastic  flow  changes  at 
points  within  the  plastic  cone.  We  have  mentioned  earlier  stepwise  linear  formu¬ 
lations  by  Pope1  and  Yamada  et  al  .2  which  employ  a  stiffness  based  upon  the  struc¬ 
tural  solution  at  the  beginning  of  the  step,  and  restrict  load  changes  so  that  the 
evolution  of  the  elastic-plastic  boundary  is  properly  traced.  In  these  formula¬ 
tions  there  is  no  account  taken  of  the  changing  constitutive  relationships  within 
the  plastic  zone  during  the  step.  Marcal  and  King3  devised  an  average  stiffness 
formulation  which  alleviates  the  need  to  regulate  load  according  to  elastic-plastic 
boundary'  movement.  It  was  designed  to  allow  an  arbitrary  step  size.  Yielding 
within  the  step  is  accommodated  by  an  element  stiffness  generated  from  an  average 
constitutive  matrix  having  weighted  elastic  and  elastic-plastic  factors.  As  in 
the  other  formulations  no  account  is  taken  of  changing  constitutive  relationships 
within  the  plastic  zone  and,  as  a  result,  the  solution  may  violate  the  yield  cri¬ 
terion  and  the  strain-hardening  law. 

Rice  and  Tracey4  considered  the  restrictions  imposed  by  the  yield  criterion 
on  stress  changes  for  a  load  step.  They  concentrated  on  the  nonhardening  model 
and  proposed  an  average  flow  rule,  corresponding  to  a  secant  approximation  to  the 
Mises  yield  surface,  which  should  be  used  in  defining  the  stiffness  so  that  the 
solution  will  not  violate  the  yield  criterion.  In  the  context  of  the  Marcal  and 
King  formulation,  implementation  requires  the  use  of  an  average  matrix  for  the 
elastic-plastic  portion  of  the  stiffness  matrix.  For  isotropic  strain-hardening 
materials,  Tracey5  extended  these  considerations  further  by  employing  step  average 
hardening  rates  for  satisfaction  of  a  given  hardening  law. 

In  each  of  these  latter  formulations3"5  the  stiffness  is  defined  in  terms  of 
quantities  which  depend  upon  the  nodal  solution.  In  the  notation  of  the  refer¬ 
ences,  the  scalar  m  for  weighting  the  elastic  and  elastic-plastic  factors  of  the 
stiffness  matrix  depends  upon  AU^;  the  same  is  true  for  the  average  flow  vector 
n,  and  the  average  hardening  rate  Y'.  Hence  these  formulations  are  stepwise 
nonlinear  and  have  an  equilibrium  equation  which  takes  the  form  of  (3),  Although 
the  averaging  methods  provide  a  solution  consistent  with  basic  constitutive  re¬ 
quirements,  the  solution  nonetheless  is  approximate.  In  fact,  only  if  the  actual 
stiffness  is  constant  for  the  step  would  the  solution  be  exact.  We  next  consider 
the  issues  involved  in  regulating  the  level  of  approximation  for  the  nonhardening 

4.  R1CF.,  J.  R.,  and  TRACEY,  D.  M.  Computational  Fracture  Mechanics  in  Numerical  and  Computer  Methods  in  Structural 
Mechanics,  S.  J.  Fenves  et  al.,  ed..  Academic  Press,  New  York,  1973,  p.  585-623. 

5.  TRACEY,  D.  M.  On  the  Fracture  Mechanics  Analysis  of  Clastic-Plastic  Materials  Using  the  Finite  Clement  Method.  Ph.P.  Thesis, 
Brown  University,  Providence,  Rhode  Island,  1973. 
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formulation  which  employs  both  the  Marcal  and  King3  arbitrary  load  change  stiff¬ 
ness  averaging,  and  the  average  flow  rule  for  yielded  points:  the  formulation 
discussed  by  Rice  and  Tracey. 

When  discussing  the  stiffness  approximations,  it  appears  that  the  most  impor¬ 
tant  issue  to  consider  is  the  changing  elastic-plastic  boundary  during  an  arbitrary 
load  interval.  Of  course,  since  the  structural  stiffness  is  formed  from  local 
stiffnesses,  this  issue  is  treated  relative  to  the  behavior  of  distinct  material 
points.  In  general,  a  point  can  experience  a  number  of  different  stages  of  de¬ 
formation  during  an  arbitrary  load  interval,  corresponding  to  elastic  response 
to  yield,  followed  by  various  phases  of  plastic  deformation,  elastic  unloading, 
and  reyielding.  The  simplest  interval  history  would  consist  of  solely  elastic 
behavior  with  the  stress  state  always  below  yield.  In  this  case  the  local  stiff¬ 
ness  remains  constant.  A  more  involved  history  is  depicted  in  igure  1.  The  de- 
viatoric  stress  space  plot  illustrates  the  important  states  for  i  point  which 
begins  a  load  interval  below  yield  at  stress  deforms  elastically  until  the 
incipient  yield  state  Sj,  and  then  undergoes  continuous  elastic-plastic  deforma¬ 
tion,  reaching  a  final  yield  surface  stress  state  Sf.  For  this  case  the  local 
stiffness  changes  abruptly  when  Sf  is  reached,  and  thereafter  it  changes  gradually 
as  the  stress  state  follows  the  yield  surface  to  Sf.  A  complex  history  entailing 
more  than  a  single  stage  of  plastic  deformation  with  elastic  unloadings  would  be 
represented  in  the  plot  by  distinct  stress  excursions  along  and  inside  the  yield 
surface  from  Sj  to  Sf. 

The  Marcal  and  King  average  stiffness  procedure  cannot  accurately  accomodate 
a  complex  interval  history.  In  fact,  the  procedure  is  strictly  correct  (m  and  Sj 
are  exact)  only  when  the  displacement  rate  of  change  is  constant  throughout  the 
load  interval.  The  solution  history  from  ^  to  Sf  is  being  approximated  by  All  ■ 
whenever  the  displacement  rate  is  nonuniform.  Hence  for  complex  histories  it  is 
clear  that  totally  unsatisfactory  approximations  can  result  through  the  implicit 
neglect  of  elastic-plastic  phases  of  the  load  interval. 
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From  the  above  it  is  clear  that  in  incremental  elastic-plastic  analysis,  load 
steps  should  be  restricted  to  encompass  portions  of  the  structural  solution  history 
which  involve  mildly  varying  displacement  rates  in  plast ical lv  deforming  regions. 
Beyond  this  there  is  the  further  need  to  restrict  load  so  that  the  average  flow 
rule  reasonably  represents  the  plastic  deformation  history.  This  can  be  accom¬ 
plished  by  restricting  the  stress  change  S^-Sj  which  occurs  during  plastic  defor¬ 
mation.  Yet,  perhaps  more  significant  is  the  fact  that  by  restricting  this  stress 
change,  displacement  rate  nonuniformities  present  within  the  load  step  are  also 
regulated.  Hence,  we  have  identified  a  single  field  variable  which  directly 
governs  the  level  of  stiffness  approximation  for  our  problems. 

Our  constraint  condition  follows  from  the  above  considerations,  and  it  in¬ 
volves  the  modulus  of  S^.-Sj,  denoted  as  ASsec  below.  At  each  load  step  we  seek 
that  solution  which  has  a  maximum  ,\Sscc  value  equal  to  a  specified  fraction  of 
the  yield  stress  Y.  If  a  represents  the  specified  fraction  then  111  takes  the 
form 

g  (AUj)  =  -  uY  =  0.  (10) 

We  have  not  attempted  to  establish  the  relationship  between  structural  stiffness 
approximation  level  and  the  constraint  parameter  a.  For  the  general  problem  this 
does  not  appear  to  be  possible.  Krieg  and  krieg''  have  considered  a  related  ques¬ 
tion.  They  discuss  yield  surface  stress  error  for  homogeneous,  constant  displace¬ 
ment  rate  deformation  using  several  constitutive  approximations  including  the 
average  normal.  We  can  conclude  only  that  convergence  to  the  exact  solution  can 
be  achieved  with  decreasing  a  values.  The  test  examples  of  the  next  section 
demonstrate  the  approach. 

NUMERICAL  RESULTS 

Here  we  discuss  two  elastic-plastic  problems  which  were  solved  using  the 
variable  load  step  approach.  The  pertinent  geometric  and  loading  features  of  the 
problems  are  shown  in  Figure  2.  Following  the  previous  development  the  nonharden¬ 
ing  Prandt 1 -Reuss  constitutive  idealization  was  used,  as  was  ti:e  solution  con¬ 
straint  condition  (10)  in  conjunction  with  the  tangent  modulus  formulation 
described  by  Rice  and  Tracey.  In  the  discussion  below,  Y  represents  the  yield 
stress  in  simple  tension  and  F.  represents  the  elastic  Young's  modulus.  Solutions 
were  obtained  using  a  Poisson's  ratio  equal  to  0.3. 

First  we  consider  a  plate  in  plane  strain  which  is  under  imposed  uniaxial 
extension.  Hill6 7  has  given  the  exact  solution  to  this  biaxial  stress  problem. 

An  interesting  aspect  of  the  solution  is  the  load-extension  relationship  from 
incipient  yield  to  limit  load.  Whereas  only  a  slight  increase  in  applied  tension 
is  possible  after  yielding  and  before  uncontrolled  plastic  deformation  occurs, 
the  displacement  necessary  to  reach  this  limit  load  state  is  unbounded.  Hence 


6.  KRIEG,  R.  D.,  and  KRIEG.  D,  B.  Accuracies  of  Sumerical  Solution  Methods  for  the  l-'lastic-Perfectlv  Nasiic  Model.  Trans. 
ASME,  J.  Pressure  Vessel  Tech..  November  1977,  p.  510-515. 

7.  HILL,  R.  The  Mathematical  Theory  of  Plasticity.  Claredon  Press.  Oxford,  England,  1950.  p.  77-"9. 
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a.  Plane  strain  uniaxial  extension  elastic-plastic  test  problem. 
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b.  Plane  stress  uniaxial  tension  of  elastic-plastic 
square  sheet  with  circular  cutout  test  problem. 

Figure  2. 


this  provides  a  valuable  test  case  for  our  solution  approach.  There  are  no  spa¬ 
tial  variations  to  contend  with  in  the  problem  so  that  a  single  low-order  element 
adequately  models  the  plate.  A  four-node  isoparametric  element  was  used  in  the 
analysis. 

In  Figure  3  four  numerical  solutions  are  given  along  with  the  exact  solution. 
The  solutions  correspond  to  a  values  of  0.2,  0.1,  0.05,  and  0.025,  where  a  is  the 
freely  specified  constraint  parameter  in  (10).  The  results  are  presented  in  plots 
depicting  imposed  stress/yield  stress  versus  imposed  strain/yield  strain,  (P/T)/Y 
versus  (U/L)/(Y/E).  These  solutions  result  regardless  of  whether  P  or  U  is  taken 
as  the  independent  loading  parameter.  The  numerical  data,  labeled  with  their 
associated  step  numbers,  are  connected  to  form  piecewise  linear  approximations  to 
the  exact  solution  which  is  given  by  the  dashed  curve  on  each  plot.  As  would  be 
expected,  the  approximations  improve,  the  step  sizes  decrease,  and  the  number  of 
steps  to  final  load  increases  as  a  is  reduced.  All  solutions  were  generated  by 
specifying  a  convergence  tolerance  value  6  equal  to  10"5.  Three  cycles  were  re¬ 
quired  to  meet  this  convergence  test  at  each  step  of  each  solution. 

A  limit  load  detection  test  was  employed  after  each  iteration  cycle  in  the 
load  step  solution  algorithm.  If  a  material  point  was  found  which  had  an  equiva¬ 
lent  plastic  strain  increase  more  than  one  thousand  times  greater  than  its  devia- 
toric  elastic  strain  modulus  change,  execution  was  terminated.  Such  a  situation 
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a  ~  005  a  =  0.025 


Rgure  3.  Load-extension  results  from  incipient  yield  to  limit  load  for 
plane  strain  example  for  a  values  from  0.2  to  0.025. 

would  approximate  the  deformation  at  limit  load  which  is  purely  plastic  and  di¬ 
rected  normal  to  the  yield  surface  at  the  limiting  stress  state.  As  indicated  in 
Figure  3,  this  test  was  met  at  load  steps  3,  4,  6,  and  9  for  the  decreasing  a 
values.  Furthermore,  the  test  was  met  in  each  case  at  an  applied  tension  very 
close  to  the  exact  limiting  value  of  1.155Y. 

The  actual  number  of  steps  required  for  each  a  choice  is  predictable  for  this 
problem.  The  constraint  parameter  restricts  deviatoric  stress  change  during  plas¬ 
tic  deformation  to  a  value  equal  to  aY  and  for  this  case  the  actual  change  is 
known  to  equal  0.185Y.  This  corresponds  to  a  biaxial  stress  change  from  (1.125Y, 
0.338Y)  at  incipient  yield  to  (1.155Y,  0.577Y)  at  limit  load. 

This  problem  served  to  elucidate  an  important  aspect  of  numerical  solutions 
near  limit  load  which  are  obtained  using  the  tangent  modulus  approach.  By  virtue 
of  the  chosen  constraint  parameter  a  (or  load  step  in  the  fixed  load  approach)  a 


solution  can  be  sought  and  found  which  places  the  stress  state  beyond  its  exact 
limiting  state.  For  the  present  problem  this  was  observed  at  step  2  of  the 
a  =  0.2  solution.  The  resulting  biaxial  stress  state  (1.1545Y,  0.S959Y)  satisfies 
the  yield  criterion,  but  it  represents  a  point  on  the  yield  surface  beyond  the 
limit  state.  This  type  of  solution  is  achievable  because  an  approximate  stiff¬ 
ness  is  being  employed.  Nonetheless,  there  is  no  insurmountable  difficulty  with 
this  behavior  since  the  error  can  be  reduced  by  decreasing  a  (or  step  size). 

The  plot  which  dramatically  illustrates  the  worth  of  the  variable  load  step 
approach  is  given  in  Figure  4.  The  plot  gives  the  a  =  0.025  discretization  re¬ 
sults,  in  the  form  of  step  size  versus  step  number,  for  both  the  force  P  and  the 
displacement  U  loading  conditions.  Unless  one  has  a  detailed  knowledge  of  the 
exact  solution,  the  discretizations  have  an  unexpected  character,  suggesting  that 
this  problem  would  entail  involved  trial  and  error  reanalysis  with  the  standard 
fixed  load  approach.  The  results  are  displayed  relative  to  Pj  and  Uj,  the  applied 
force  and  displacement  at  incipient  yield.  The  force  boundary  condition  AP/Pj 
data  (marked  with  triangles)  is  plotted  using  the  left  axis  scale.  The  right  axis 
applies  to  the  displacement  boundary  condition  AU/Uj  data  (marked  with  X's).  The 
AP/Pj  values  vary  from  0.0066  for  step  2,  the  first  step  after  incipient  yield  is 
reached,  to  0.0009  at  limit  load.  The  corresponding  AU/Uj  values  range  from 
0.136  to  0.940. 


Figure  4.  Load  discretization  results  for  plane  strain  problem  with 
constraint  parameter  a  equal  to  0.025. 
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The  second  problem  involves  the  plane  stress  uniaxial  tension  of  a  square 
sheet  with  a  centered  circular  cutout.  This  problem  is  more  typical  of  those  en¬ 
countered  in  practical  analysis  in  the  sense  that  there  are  undoubtedly  complicated 
spatial  variations  in  the  structure  which  significantly  change  character  as  the 
yielding  progresses.  However,  a  comprehensive  spatial  convergence  study  was  not 
undertaken,  for  as  throughout  this  paper  the  emphasis  was  placed  on  load  discretion 
tion  and  how  it  affects  the  solution  for  an  arbitrary  finite  element  model.  Two 
models  consisting  of  four-node  isoparametric  elements  were  employed:  one  had  a 
3x6  element  arrangement  (28  nodes  1,  while  the  other  had  a  5x8  (54  node)  arrangement 

Load-extension  results  are  given  in  Figure  5  for  «  choices  of  0.05  and  0.15. 
The  data  are  plotted  in  the  normalised  form  (P/H)/Y  versus  (V/H)/ (Y/F.) ,  where  P/H 
is  the  uniform  tension  across  the  ends  of  the  plate,  and  V  is  the  displacement  of 
the  center  of  the  loaded  edge.  The  step  data  are  numbered  for  the  a  =  0.15  solu¬ 
tion  and  the  spread  of  the  plastic  sone  is  illustrated  by  the  shaded  elements  in 
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Figure  5.  Load-extension  results  from  load  free  state  to  limit  load  for 
plane  stress  example  for  a  values  of  0.15  and  0.05. 


the  sketches  of  the  model.  As  can  be  seen,  the  significant  differences  in  the 
two  solutions  begin  at  the  knee  of  the  curve  when  yielding  loses  its  localized 
character. 


Limit  load  was  not  detected  numerically  by  the  above-mentioned  strain  test, 
hence  loading  was  continued  until  step  size  was  reduced  to  very  small  fractions 
of  the  incipient  yield  load  Pj.  The  discretization  results  are  plotted  in 
Figure  6.  As  in  the  previous  problem  it  is  unlikely  that  a  priori  judgment  would 
suggest  the  form  of  the  results,  with  L?/?i  starting  at  values  of  0.447  and  0.188 
and  ending  at  values  in  the  neighborhood  of  0.009  and  0.001  for  o  equal  to  0.15 
and  0.05,  respectively. 

A  limit  analysis  of  this  problem  has  been  presented  in  the  literature.8  It 
places  the  exact  limit  value  of  (P/H)/Y  between  0.31  and  0.40,  whereas  our  limit 
load  result  for  a  =  0.15  predicts  a  value  of  0.46.  A  solution  using  the  5x8  mesh 
with  a  =  0.15  was  obtained,  and  it  predicts  a  limit  load  value  of  0.44.  This  in¬ 
dicates  that  spatial  discretization  is  the  source  of  the  observed  differences. 


This  problem  provided  useful  data  concerning  the  convergence  properties  of 
t  e  solution  algorithm  at  a  load  step.  As  would  be  expected,  it  was  found  that 
t  e  typical  number  of  cycles  to  meet  the  6  test  increases  with  a  value,  and  fur- 
thermore  the  required  number  increases  at  steps  near  limit  load.  A  6  value  of 
10  was  used  for  this  problem.  Cycle  counts  averaged  close  to  5  for  a  =  0  05 
ana  in  the  neighborhood  of  10  for  a  =  0.15.  ' 

As  mentioned  earlier,  the  solutions  were  obtained  using  a  formulation  which 
employs  an  average  yield  surface  normal4  to  define  the  plastic  flow  rule  for 
yielded  points.  Special  considerations  were  necessary  in  adapting  this  averaging 
technique  to  the  plane  stress  problem.  In  plane  stress  there  is  the  need  to 
establish  the  average  normal  in  terms  of  the  out-of-plane  direct  strain  increment, 
but  this  strain  component  depends  upon  the  flow  rule  for  its  definition.  A  method 
was  devised  for  defining  these  quantities  in  a  way  which  insures  that  both  the 

P*a"aJ  c°ndlt^n  aT}d  the  yield  condition  are  satisfied.  The  details  of  the 

method  will  be  described  in  a  forthcoming  report. 

CONCLUSIONS 

WhiCh  n°W  include  the  load  Path  discretization  and 
corresponding  field  solution,  demonstrate  the  viability  of  our  variable  load  step 

solution  algorithm  in  elastic-plastic  analysis.  We  expect  that  the  algorithm  will 
apply  equally  well  to  other  nonlinear  problems  treated  by  stepwise  nonlinear  tan- 

f^Ll?0dUlUS-50rmaia^10nS’r- IC  is  clear  that  the  success  of  the  approach  is  pred¬ 
icated  upon  identifying  a  field  variable  which  controls  the  level  of  stiffness 

and  s^itab;y  constraining  the  variable  at  each  step  of  the  solution. 
"  1  pf  the  unPredictable  character  of  the  results  obtained,  the  customary  fixed 

J°*d t  tl  approach  now  appears  tenuous.  The  variable  step  approach  not  only  elim¬ 
inates  the  requirement  for  a  pnori  discretization  but  also  provides  a  series  of 
consistent  incremental  solutions  according  to  the  desired  level  of  approximation 
The  analyst  need  only  specify  the  proportional  loading  segments  of  the  load  path 
and  supply  values  for  the  constraint  parameter  (a  in  our  elastic-plastic  problems) 
a"d  tb*  Convergence  tolerance  parameter  6  for  the  iterative  solution  at  each  step. 
With  these  parameters  the  analyst  can  efficiently  examine  and  control  solution 

equilibrium^  ^  regardS  stlffness  approximation  level  and  overall  satisfaction  of 
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